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A fast simulation algorithm for the calculation of multitime correlation functions of open quan- 
tum systems is presented. It is demonstrated that any stochastic process which "unravels" the 
quantum Master equation can be used for the calculation of matrix elements of reduced Heisen- 
berg picture operators, and thus for the calculation of multitime correlation functions, by extending 
the stochastic process to a doubled Hilbert space. The numerical performance of the stochastic 
simulation algorithm is investigated by means of a standard example. 
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I. INTRODUCTION 

The state of an open quantum system is conventionally 
described through a reduced density matrix p(t) whose 
dynamics is given by a dissipative equation of motion - 
the quantum Master equation. From a numerical point of 
view this formalism has a major drawback: for a system 
whose state is described in a TV-dimensional Hilbert space 
H, the quantum Master equation is a set of N(N + I )/2 
coupled differential equations. Hence the numerical eval- 
uation of the quantum Master equation is in practice not 
feasible for large systems [Q. 

This difficulty does not arise in the stochastic wave 
function approach to open systems 0-0]: Here, the state 
of an open quantum system is described by a stochas- 
tic wave function ijj{t) £ Ti., i.e., by a TV— dimensional 
state vector. The stochastic time evolution of ip(t) is ei- 
ther defined through a stochastic Schrodinger equation 
§,§ (which is a stochastic differential equation) or al- 
ternatively through a conditional transition probability 
T[ip, t\ipo, to] 1^,0, which is the probability density of 
finding the system in the state ifi at time t > to under the 
condition that the system is in the state ipo at time to- 
The connection to the density matrix formalism is made 
through the relation 



(1) 



p(t) = J DybDr J DrboDro MM 
xT[ip,t\il> ,t ]P[ip ,t }, 



where P["0Oj^o] is the probability density of an ensem- 
ble of normalized pure states characterizing some initial 
density matrix po and DipDip* is the Hilbert space vol- 
ume element ||[rj . The integrals extend over the Hilbert 
space Ti. This relation ensures that one-time expectation 
values of any system operators are calculated correctly. 
Note that this condition alone does not uniquely specify 
a stochastic process. Diffusion type stochastic processes 
[B,p| as well as piecewise deterministic jump processes 
\2-yj have been proposed in the literature. A unique 
stochastic process can only be derived by making further 



assumptions such as specifying a certain measurement 
scheme |^||. 

Especially in quantum optical systems one-time expec- 
tation values of system observables are not the only mea- 
surable quantities: for example, the spectrum of fluores- 
cence of a two level system is the Fourier transform of the 
two-time correlation function ((ct + (t)<7~)) s in the station- 
ary state, where cr ± denote the pseudo spin operators of 
the system. Thus, for a complete description of open 
quantum systems it is necessary to introduce Heisenberg 
picture operators. In the density matrix formalism this 
concept is well understood [|k|: Consider the quantum 
Master equation 



p(t) = L(t)p{t), 



(2) 



where the super-operator L(t) is defined as 

L(t)p(t) = -i[H(t) >P (t)] (3) 

i 

The operator H(t) is essentially the Hamiltonian of the 
isolated system which contributes to the coherent part of 
the dynamics, and the rates "fi and the Lindblad oper- 
ators Ji describe the dissipative coupling of the system 
to its environment through the i-th decay channel. The 
solution of eq. (||) with respect to the initial condition 
p(to) — po can be expressed for t > to in terms of the 
propagation super-operator V(t, to) as [ flij 



p(t) = V(t,t )po, 



(4) 



where V(t, to) is the solution of the differential equation 
d 



dt 



V(t,t ) = L(t)V(t,t ), 



(5) 



with the initial condition V(to,to) = I. For an arbitrary 
Schrodinger system operator A the matrix elements of 
the reduced Heisenberg picture operator are defined as 



i 



A t (0 o ,Vo) = (^o,*o|^(*)|^o,*o> 



= Tr 



ays 



{AV(t,to)\ifo)(<l>o\}. (6) 



Eq. (Q) can be interpreted in the following way: for the 
calculation of the matrix element A t ((j} ,ijj ) start with 
the initial "density matrix" |'i/'o)( < ^o| an( i propagate it up 
to the time t. Then calculate the expectation value of A 
with respect to the propagated "density matrix" . How- 
ever, since |^o) (<^o I is in general not a positive matrix and 
thus not a true density matrix, it can not be character- 
ized by a probability density P[ipo, to] of normalized pure 
states in Ji (cf . eq. (|l]) . Hence a direct application of the 
stochastic wave function approach to the calculation of 
Heisenberg picture operators is not possible. 



II. HEISENBERG PICTURE OPERATORS IN 
THE STOCHASTIC WAVE FUNCTION 
APPROACH 



For the piecewise deterministic jump process the sim- 
ulation algorithm reads as follows: 1) Start with the nor- 
malized state 9q — (0o, V'o) T /v / 2 at to. 2) Draw a random 
number r}\ from a uniform distribution on [0, 1]; this ran- 
dom number will determine the time of the first jump. 3) 
Propagate 9q according to the Schrodinger-type equation 



i—6{s) = H cB (s)9(s) 
as 



(8) 



where the extensions of the Hamiltonian H and Lindblad 
operators Ji to the doubled Hilbert space are defined as 



H(s) = 



_ ( H{s) 



H(s) 



Ji = 



J l 
Ji 



(9) 



and the non-Hcrmitian effective Hamiltonian is defined 



H eS {s)=H{s)- l -Y,^ J i J i- 



(10) 



In a closed system where the time evolution of states is 
given through the unitary propagator U(t, to) we can cal- 
culate arbitrary matrix elements A t {4>o, ipo) of a Heisen- 
berg operator A(t) in the following way (cf. fig. [I]): 
propagate <f>o and ipo to obtain <fi — U(t,to)(f>o and 
ip = U(t, to)ipo, respectively and then evaluate the scalar 
product (0|j4|i/'). This method is easily generalized to 
the calculation of matrix elements of a reduced Heisen- 
berg picture operator, i.e., to open systems: Instead of 
propagating the state vectors <j>o € H and ipo € "H sepa- 
rately we can construct a stochastic process in the dou- 
bled Hilbert space TL = TL TL which propagates the 
normalized pair of state vectors 8o — {4>o, V'o) T /v / 2 G TL 
simultaneously in such a way that the following condition 
holds: 



A t (<fo,ipo) 



D9D6*{<t>\A\ij>)T[6,t\6o,to] 



(J) 



where 9 — (</>, i/') T and we introduced the conditional 
transition probability T for the stochastic process in the 
doubled Hilbert space TL. Throughout this letter, the 
superscript T denotes the transpose of a vector. This 
condition states that matrix elements of arbitrary Heisen- 
berg operators are calculated correctly. It is important 
to note that eq. (0) alone does not specify the stochastic 
time evolution in the doubled Hilbert space uniquely. In 
fact, each stochastic process which can be used to simu- 
late the quantum Master equation (^) can be extended 
to the doubled Hilbert space and used for the calculation 
of the matrix elements of arbitrary Heisenberg picture 
operators. We will first demonstrate this for the piece- 
wise deterministic jump process proposed in [§-[7| and 
then generalize the result. A derivation of the simula- 
tion algorithm for the stochastic time evolution in the 
doubled Hilbert space which is based on a microscopic 
system-reservoir model can be found in Ref. |15| . 



4) The time T\ of the first jump is determined by the 
condition 



m = \\o(Ti)\\ 



(ii) 



5) Select a particular type of jump with probability 
liWij J2ili w i> where Wi = \\ J s ;#(Ti)|| 2 . The state of the 
system immediately after the first jump is given by 



9(T 1 ) = J l 9(T 1 )/\\J l 9(T 1 )\\ 



(12) 



6) Draw a second random number 772 to determine the 
time T2 of the next jump and propagate 9{T\) according 
to the differential equation (||) and so on until s = t. 7) 
The state of the system at time t is given by 



m = m),m) T =o{t)/\\Kt)\\ 



(13) 



The matrix elements of the reduced Heisenberg picture 
operator A(t) are then obtained by computing 



(14) 



where the angular brackets ((• • •)) denote the average over 
the realizations of the stochastic process. 

In order to show that this algorithm leads to the cor- 
rect result, we introduce the density matrix 



m = 



Pu{t) P12CO 



(15) 



on the doubled Hilbert space TL which is a solution of the 
extended quantum Master equation 

?(t) = -i [£(*),?(*)] (16) 
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with the initial condition 

By definition of the extended operators H and Ji each 
component pij(t) which is an operator on TL is a solu- 
tion of the original quantum Master equation (Q) . Since 
7>2i{t) is a solution of eq. (|^) with the initial condition 
P2i(*o) = l^o)(^o|/2 the matrix elements A t (</> , ?/>o) of a 
reduced Heisenberg picture operator A(t) can be written 
as (cf. eq. (^)) 

A t (0 o ,Vo) = 2Tr sys {Ap 21 (t)}. (18) 

Now consider a particular "unraveling" of the extended 
quantum Master equation ( |l6| ) which is characterized by 
a conditional transition probability T[9,t\9o,t ] in the 
doubled Hilbert space. For the density matrix p(t) we 
then obtain in analogy to eq. (|Tj) 

p(t) = J D9D9* \9)(8\ T[6,t\e ,to]: ( 19 ) 

and hence for p2i(t) 

hi{t) = J D9D0* \i>)(4,\ T[6,t\e ,t ], (20) 

where 9 = (6, ip) T ■ By inserting eq. (|Tj|) into eq. ( |l8| ) we 
recover eq. QTj) . Thus we have shown that matrix elements 
of reduced Heisenberg picture operators are calculated 
correctly (i.e., eq. (0) holds) if the stochastic process in 
the doubled Hilbert space can be used to simulate the 
extended quantum Master equation ([l6|). Since this is 
the case for the simulation algorithm presented above we 
have completed the proof. 

It is important to note that the above proof does not 
rely on a specific "unraveling" of the quantum Master 
equation (|l6|) . On the contrary, it is valid for any stochas- 
tic process the covariance matrix (see eq. (fiT^)) of which 
is governed by eq. (|l6|). 

III. MULTITIME CORRELATION FUNCTIONS 

The simulation algorithm in the doubled Hilbert space 
can also be used for the calculation of multitime cor- 
relation functions. Consider for example the two-time 
correlation function 

g(<t>oMM) = (fo,fo|A(t2)B(ti)|0o,*o), (21) 

where t\ <t 2 - Here, the stochastic simulation algorithm 
would read as follows: 1) Start in the state <\>o at time 
to and use the stochastic time evolution in the Hilbert 
space TL to obtain the stochastic wave function 4>(ti). 2) 
Propagate the state 



9(h) = (0(i 1 ) J fl#t 1 )) T /||(0(t 1 ) J fl#t 1 ))|| (22) 

using the stochastic time evolution in the doubled Hilbert 
space to obtain the state vector 9(t 2 ) — ((j>(t 2 ), ip(t 2 )) T - 
The multitime correlation function is then obtained by 
computing 

ff(to,*l,*2) = ((||(^l), Wl))|| 2 <^2)|A|^(i 2 )))). 

(23) 

The generalization of this scheme to the calculation of 
arbitrary time-ordered multitime correlation functions of 
the form 

g(4>o,t (s ,t 1 , ...,i„,si, ...,s m ) (24) 
= (0Oj*o|4i(*i) • • ■ A n (t n )B m (s m ) ■ ■ ■ Bi(s 1 )\(f> ,t }, 

where to < ■ ■ ■ < t n , and to < si < • • • < s m , and Ai 
and Bi are arbitrary system operators is straightforward: 
Order the set of times {t\, ■ ■ ■ t n , s\, ■ ■ ■ s m } and rename 
them Ti such that n < • • • < r q where q is the number 
of distinct time points. Then define a set of Schrodinger 
operators FJ and Gi as 

!Fi = A\, Gi — I, if n = ti 7^ Sj for some i and all j, 
Gi = I,Fi = Bj, if n = sj =/= ti for some j and all i, 
Gi = A';,Fi = Bj, if n = ti = Sj for some i and j. 

(25) 

The multitime correlation function g(4>o, to, t\, 
t n , s\, s m ) is then obtained in the following way: 

1) Start with the state 4>o a t time to and propagate it 
up to the time r\ to obtain (j)(r{). 

2) Propagate the state 

9( n ) = (F 1 0(r 1 ),G 1 0(r 1 )) T /||(F 1 0(r 1 ),G 1 0(r 1 ))|| (26) 

to obtain 9(r 2 ) = (^(ra), ip(r2)) T ■ 

3) Jump to the state 

9(r 2 ) = (F 2 0(r 2 ),G 2 ^(r 2 )) T /||(F 2 0(r 2 ),G 2 V(r 2 ))|| (27) 

and propagate it up to and so on. 

g(cj>o,to,h, ...,t n , si, ...,s m ) is then given by 

g(4>o,to,ti, tn, si, s m ) = (28) 
((|| (FnKn), GnKn)) \\ 2 \\ (Fa^ra), G 2 ^(r 2 )) f • • • 

x||(F g „ 1 0(r g „ 1 ),G g „ 1 ^(r 9 _ 1 ))|| 2 <^(r 9 )|F 9 tG 9 |V((r 9 )))). 

It is important to note, that also for higher order cor- 
relation functions, we only have to propagate two state 
vectors. 

Finally, let us remark that the choice of the initial con- 
dition ( p3[ ) (or (pj(f), respectively) is not unique. We can 
also multiply the operator B by a constant e and define 
the state vector 9 E (ti) as 
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6 £ {ti) = (^(ti),eWi)) T /||(0(ii),eWi))|| (29) 
and accordingly the correlation function g((f>Q,ti,t2) as 

g(<t>o,h,h) = (30) 

' // '|(^ 1 ),ei?0(i 1 ))|| 2 (^(t 2 )|A|^(i2)) 



where 9 e (t2) — (4 > e(t2),tpe(t2)) T is obtained by propa- 
gating 9 e (t\) according to the simulation algorithm in 
the doubled Hilbert space. Again, the unnormalized de- 
terministic motion is governed by the equation of motion 



i-<t>e{t) = H BS {t)<t> e {t) 
at 

at 

but in the limit e — > 0, we find 

\Mt)\\ -> \\l(t)\\ 

w t = || lA(T)\\ 2 - \\Ji4> e {T)\\ 2 , 



(31) 
(32) 



(33) 
(34) 



and hence the jumps of the trajectory 9 E (t) are com- 
pletely governed by the jumps of 4> e (t), which evolves ac- 
cording to the "usual" stochastic time evolution in Ti (cf. 
Eqs. and @). In this limit we obtain a procedure 
first proposed by Dum et al. in Ref. which is based 
on "probing the system with 8 kicks" (see Appendix D of 
Ref. Q). For further discussions of this method see for 
example the Refs. Jl6Ml8||. 



IV. NUMERICAL RESULTS 

In order to investigate the numerical performance of 
our simulation algorithm, we compare it with the method 
proposed by Dum et al. in Ref. S and with an alterna- 
tive method proposed by Dalibard et al. which is based 
on a decomposition of the stochastic trajectory into four 
sub-trajectories 0. Note that all procedures are fully 
consistent with the quantum regression theorem [^Jl9| 
and hence lead to the same result for the multitime cor- 
relation function. However, the numerical performance 
of the algorithms is quite different. We demonstrate this 
by means of a standard example of quantum optics - the 
calculation of the spectrum of resonance fluorescence of 
a two level system. In fig. |^ (a) - (c) we show the com- 
putational time necessary to achieve a given accuracy 
(measured by the relative error of the correlation func- 
tion ((<7 + (t)<j~)) s in the stationary state) for a coherently 
driven two level atom with Rabi frequency £1 = IO7 ob- 
tained on a RS6000 workstation. The solid lines repre- 
sent the mean square deviation of the numerical solution 
from the exact solution |^| and the dashed lines show 
the mean estimated standard deviation of the numerical 
solution. Obviously, the latter quantity provides for all 
algorithms a very good measure of the accuracy of the 



numerical simulation. In fig. (d) we compare the esti- 
mated standard deviation for the three algorithms. Ob- 
viously, the numerical performance of the algorithms pro- 
posed by Dum et. al. and our algorithm is quite similar, 
although the convergence of our algorithm is smoother. 
On the other hand, for a given accuracy the stochastic 
simulation in the doubled Hilbert space is by a factor of 3 
faster than the algorithm proposed in [|| . We expect this 
result to be even better for higher order correlation func- 
tions since for a multitime correlation function of the type 
of eq. ( p4[ ) one has to propagate in general 4™+ m - 1 differ- 
ent state vectors in each realization using the method of 
Castin et. al, whereas in our approach it is only necessary 
to propagate two state vectors. 

Let us briefly summarize the main results of this letter: 
We have shown that starting from a stochastic simulation 
algorithm for the quantum Master equation (^) it is pos- 
sible to obtain a fast simulation algorithm for the calcu- 
lation of matrix elements of arbitrary Heisenberg picture 
operators and time-ordered multitime correlation func- 
tions by making the substitutions 



i/j e n — > e e h, h — ► h, ^ 



J, 



(35) 



i.e., we replace the stochastic wave function tp(t) by a 
stochastic wave function 9(t) in the doubled Hilbert space 
and extend accordingly the operators H and Ji which are 
present in the quantum Master equation to the doubled 
Hilbert space (cf. eq. (g)). We emphasize that these re- 
placements can be done for any "unraveling" of the quan- 
tum Master equation, e.g., also for the quantum state 
diffusion model |||| . The resulting stochastic process in 
the doubled Hilbert space is then similar to a process first 
proposed by Gisin in However, the latter process is 
only well defined, when the initial states (j) and i^q are 
non-orthogonal, i.e., if (</>o|^o) 7^ 0. This problem does 
not occur in the ansatz presented here. 
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FIG. 2. Calculation of the first order correlation function 
((cr + (r)cr _ )) s for a coherently driven two level atom on reso- 
nance. This figure shows the relative error vs. the CPU time 
in seconds for the simulation algorithms proposed in Ref. || 
(a), Ref. H (b), and for our algorithm in the doubled Hilbert 
space (c). The solid lines represent the mean square devi- 
ation of the numerical solution from the exact solution and 
the dashed lines show the estimated standard deviation of the 
numerical solution. In Fig. ^ (d), we compare the estimated 
standard deviation for the three algorithms. 



(a) Unitary time evolution 

[/(Mo) 

► \<f>) 

(<t>\m 

m — * m 

(b) Stochastic time evolution 

T[(^),t|(^o)/V2,t ] 

\<M — 



i i i i mm) 

l^o) T7 T7 ► Wi 

FIG. 1. Calculation of Heisenberg operator matrix ele- 
ments: (a) for a closed system and (b) for an open system. 
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